09 - Remote Sensing with Landsat package

Working with satellite imagery in R

In this exercise you learn to use the functions of the landsat package to process satellite imagery, specifically, how to:

  • streamline histograms so you can combine images
  • manipulate multi-band imagery,
  • extract and create new data from images, such as NDVI and SAVI
  • classify image values using kmeans unsupervised algorithm to detect similar areas
  • segment features in satellite imagery

Task 1: Set up your workspace

Start by installing the required packages first. rasterVis and RColoBrewer are used for visualisation of rasters, lattice and latticeExtra for extra graphical utilities, while landsat provides the imagery samples and rgl allows for 3D visualisation.

library(lattice)
library(latticeExtra)
library(RColorBrewer)
library(rasterVis)
library(rgl)  # macs might need another prerequisite
library(landsat)
library(terra)

Task 2: Pre-processing of Landsat datasets

Landsat packages offers sample Landsat satellite imagery decomposed into single bands, and labelled by the month of capture (25 November 2002). You will load and practice raster analysis with these 300x300 pixel samples. Note that each single band image shows as panchromatic (greyscale) with values ranging from 0 (white) to 255 (black), with the interim colors being shades of grey. This is unsigned 8 bit imagery. The number behind the filename (e.g. nov3) refers to the color band of the image: 1 = Green, 2 = Blue, 3= Red, 4 = Near-infrared. Once you combine these color bands and plot with plotRGB(), you can see the true- or false-color imagery depending on which color you load into which RGB channel.

The images are in the SpatialGridDataFrame format and need to be converted to raster format before manipulation.

## load indvidual band image data from landsat package
`?`(nov)

# load band#3 red channel of the image
data(nov3)
plot(nov3)

data(nov4)
plot(nov4)

Task 3: Load elevation data and plot it in 3D (Adela to demo

dem in the landsat package denotes a ‘digital elevation model’. Once you load it you can convert it to Formal Raster Layer with raster() function and then continue processing using the usual raster package functions.

plot3D() is a neat function in the rasterVis package, which opens a separate windows and plots elevation values in 3D space if present in the raster. You need to close the window if you want to update or plot another object.

# Load and plot the digital elevation model in landsat package
library(raster)
data(dem)
plot(dem)

dem <- raster(dem)

Plot in 3D with raster

RasterVis package function plot3D() loads a neat 3D viewer, which opens in a new window If you have the raster library and don’t see the plot3D() function working, enable the rglwidget().

library(raster)
library(landsat)
rglwidget()
plot3D(dem, rev = T, zfac = 1)

I have not been able to get an equivalent 3D image with terra library

Task 4: Load and explore RGB data components

# let's load data for July and explore it
data("july1")
data("july2")
data("july3")
data("july4")

j1 <- rast(july1)  # blue
j2 <- rast(july2)  # green
j3 <- rast(july3)  # red
j4 <- rast(july4)  # near-infrared

## check out the image histogram\t
plot(j1)

plot(j4)

hist(j1, main = "Band 1 - Blue of Landsat")

boxplot(j1, main = "Band 1 - Blue of Landsat")

Task 5: Plot RGB image

library(terra)

# Create a true-color composite (R = j3, G = j2, B = j1)
myRGB <- c(j3, j2, j1)  # `c()` combines into a SpatRaster with multiple layers

# Create a false-color composite (CIR: NIR = j4, R = j3, G = j2)
myCIR <- c(j4, j3, j2)

### let's see how the NIR, R, and G bands relate (from lattice)
names(myCIR) <- c("NIR", "Red", "Green")
splom(myCIR, varname.cex = 1)  # scatter plot matrix

## better
plotRGB(myRGB)

plotRGB(myCIR)

Task 6: Manipulate image rendering by histogram stretching and cloud masking

First, let’s use histogram stretch

## different stretches - here histogram based
plotRGB(myRGB, stretch = "hist")

plotRGB(myCIR, stretch = "hist")

Next, let’s try linear stretch

## different stretches - here linear stretch
plotRGB(myRGB, stretch = "lin", main = "True Color RGB")

plotRGB(myCIR, stretch = "lin", main = "False Color CIR")  # in CIR red = green!

Any idea what the red color represents in myCIR?

Finally, drape one of the images over a 3D model. This bit can be a bit touchy, relies on the raster library and take a while to get to work. I tend to run the plot3D() lines alone to generate the 3D view in a pop-up window rather than rmarkdown output. Everytime you wish to refresh the view, you must close the pop-up window.

## finally...  rglwidget() # this widget helps get the first view rendered in
## rmarkdown, refreshing however is more tedious

# ?rglwidget()

t <- plot3D(dem, col = rainbow(255))  ## you need to close RGL device manually first and then run this line!


# Save it to a file.  This requires pandoc
filename <- tempfile(fileext = ".html")
htmlwidgets::saveWidget(rglwidget(), filename)
browseURL(filename)

plot3D(dem, drape = raster(july4))  ## should drape image j4 over DEM, if problematic try in .R script and watch for a pop-up widget. Ask Adela to demo!

We have a lot of clouds in the July image - mask them!

RStoolbox comes with a suite of pre-processing functions, including cloudMask to identify clouds in optical satellite imagery:

library(landsat)

j1 <- rast(july1)  # blue
j2 <- rast(july2)  # green
j3 <- rast(july3)  # red
j4 <- rast(july4)  # near-infrared

# Create a false-color composite (CIR: NIR = j4, R = j3, G = j2, B = j1)
myCIR <- c(j4, j3, j2, j1)

# Rename the NIR, R, and G bands
names(myCIR) <- c("NIR", "Red", "Green", "Blue")
# Calculate cloud index
library(RStoolbox)
library(ggplot2)

cldmsk <- cloudMask(myCIR, blue = 4, tir = 1)
ggR(cldmsk, 2, geom_raster = TRUE)  # geom_raster plots faster, 2 indicates the mask

# mask by threshold, region-growing around the core cloud pixels
cldmsk_final <- cloudMask(cldmsk, threshold = 0.1, buffer = 5)

## plot cloudmask
library(ggplot2)
ggRGB(myCIR, stretch = "lin") + ggR(cldmsk_final[[1]], ggLayer = TRUE, forceCat = TRUE,
    geom_raster = TRUE) + scale_fill_manual(values = c("yellow"), na.value = NA)

Task 7: Create new information from satellite imagery: NDVI

Let us calculate the Normalized Difference Vegetation Index (NDVI) and see where the vegetation grows most in our Landsat image: Remember the formula for NDVI is: (NIR - RED) / (NIR + RED)

# prep imagery
n3 <- rast(nov3)  # RED
n4 <- rast(nov4)  # NIR

ndvi <- (n4 - n3)/(n4 + n3)
ndvi
class       : SpatRaster 
dimensions  : 300, 300, 1  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : 390045, 399045, 4482105, 4491105  (xmin, xmax, ymin, ymax)
coord. ref. :  
source(s)   : memory
name        : /data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc 
min value   :                                               -0.3114754 
max value   :                                                0.5664336 
plot(ndvi)

# uncomment and run in 3D plot3D(dem, drape=ndvi, zfac = 1.5)

### remove values below zero
ndvi[ndvi <= 0] <- NA
plot(ndvi)

Now you can see the areas of the highest reflectance and thus most healthy vegetation in November

Task 8: Create new information from satellite imagery: SAVI

SAVI stands for Soils Adjusted Vegetation Index, and this is another calibrated view of the ground.

### another index SAVI (soil adjusted vegetation index)
ndvi <- (n4 - n3)/(n4 + n3)
savi <- (n4 - n3)/((n4 + n3) * 0.25)  # with L=1 -> similar to NDVI

### let´s compare visually
par(mfrow = c(1, 2))

plot(savi, main = "SAVI")
plot(ndvi, main = "NDVI")

par(mfrow = c(1, 1))

Task 9: Unsupervised Classification with k-means

We would like to isolate and better see the clusters of growth within our image. We will run kmeans function on a composite image in order to cluster data based on similarity or similar groups!

# first, let's select an image and make it into a brick including ndvi
data(nov2)
data(nov1)
n2 <- rast(nov2)
n1 <- rast(nov1)
ndvi
class       : SpatRaster 
dimensions  : 300, 300, 1  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : 390045, 399045, 4482105, 4491105  (xmin, xmax, ymin, ymax)
coord. ref. :  
source(s)   : memory
name        : /data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc 
min value   :                                               -0.3114754 
max value   :                                                0.5664336 
which(is.na(as.data.frame(ndvi)))  # should be 0, rerun the ndvi creation if not.
integer(0)
# create a new composite brick out of the available data
myNewBrick <- c(n4, n3, n2, n1, ndvi)
splom(myNewBrick)

plot(myNewBrick)

Next, run the kmeans classification. Beware that the kmeans function does not tolerate NA/INF/NaN and similar values. Our new brick should not have any but in future classification remember that you need to get around them, either by exclusion or substitution via mean values.

# Run kmeans classification on the values in your new brick Read on
# Thresholding here:
# https://rspatial.org/raster/rs/3-basicmath.html#vegetation-indices
ICE_df <- as.data.frame(myNewBrick)
set.seed(99)


cluster_ICE <- kmeans(ICE_df, 4)  ### kmeans, with 4 clusters
str(cluster_ICE)
List of 9
 $ cluster     : Named int [1:90000] 4 3 4 4 4 1 4 4 3 3 ...
  ..- attr(*, "names")= chr [1:90000] "1" "2" "3" "4" ...
 $ centers     : num [1:4, 1:5] 48.2 37.5 79.8 59.2 39.9 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:4] "1" "2" "3" "4"
  .. ..$ : chr [1:5] "/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc" "/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band3.geo.asc" "/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band2.geo.asc" "/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band1.geo.asc" ...
 $ totss       : num 20611557
 $ withinss    : num [1:4] 1042091 848325 839810 1215809
 $ tot.withinss: num 3946035
 $ betweenss   : num 16665522
 $ size        : int [1:4] 34968 29608 8009 17415
 $ iter        : int 3
 $ ifault      : int 0
 - attr(*, "class")= chr "kmeans"
# cluster_ICE <- cluster::clara(ICE_df, 4) ### another option, clara, with 4
# clusters

# convert cluster information into a raster for plotting
clusters <- rast(myNewBrick)  ## create an empty raster with same extent than ICE
clusters <- setValues(clusters, cluster_ICE$cluster)  # convert cluster values into raster
clusters
class       : SpatRaster 
dimensions  : 300, 300, 5  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : 390045, 399045, 4482105, 4491105  (xmin, xmax, ymin, ymax)
coord. ref. :  
source(s)   : memory
names       : /data/~eo.asc, /data/~eo.asc, /data/~eo.asc, /data/~eo.asc, /data/~eo.asc 
min values  :             1,             1,             1,             1,             1 
max values  :             4,             4,             4,             4,             4 
plot(clusters)

# uncomment to plot the clusters in 3D over the DEM plot3D(dem, drape=clusters,
# col=c('white', 'green', 'blue', 'yellow'))

# calculate the average spectral signature of 1-4 bands of growth
ICE_mean <- zonal(myNewBrick, clusters, fun = "mean")
ICE_mean  # see the values for ndvi (layer) being most distinct
  X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band4.geo.asc
1                                                         1
2                                                         2
3                                                         3
4                                                         4
  X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band3.geo.asc
1                                                         1
2                                                         2
3                                                         3
4                                                         4
  X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band2.geo.asc
1                                                         1
2                                                         2
3                                                         3
4                                                         4
  X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band1.geo.asc
1                                                         1
2                                                         2
3                                                         3
4                                                         4
  X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band4.geo.asc.1
1                                                           1
2                                                           2
3                                                           3
4                                                           4
  /data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc
1                                                 48.23805
2                                                 37.47494
3                                                 79.83131
4                                                 59.23101
  /data/LANDSAT/MASTER/Gorskiout/G20021125_7_band3.geo.asc
1                                                 39.87062
2                                                 33.39692
3                                                 41.67599
4                                                 45.38708
  /data/LANDSAT/MASTER/Gorskiout/G20021125_7_band2.geo.asc
1                                                 39.84423
2                                                 35.98109
3                                                 45.57810
4                                                 44.90479
  /data/LANDSAT/MASTER/Gorskiout/G20021125_7_band1.geo.asc
1                                                 55.70399
2                                                 53.02459
3                                                 58.27282
4                                                 58.88780
  /data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc.1
1                                                 0.09473271
2                                                 0.05560162
3                                                 0.31161705
4                                                 0.13208106

Note that you have aggregated the final cluster raster by using the focal() function using the mean function on the four clusters identified by kmeans() as similar.

Task 10: What is the trend in de-/afforestation? - Individual tree crown segmentation

The ITC (Individual Tree Crowns) delineation approach finds local maxima within imagery that contains subtle color differences, such as the canopy image provided. The itcIMG() function designates these maxima as tree tops, then uses a decision tree method to grow individual crowns around the local maxima.

The image we use is based on LiDAR (Light Detection and Ranging) in xyz format. To get the segmentation going you will need the terra and itcSegment packages

# install.packages('itcSegment')
library(terra)
library(itcSegment)

forest <- rast("../data/imgData.tif")

plot(forest)

# Use the itcIMG() function to detect and grow the individual crowns
`?`(itcIMG())
se <- itcIMG(forest, epsg = 32632)
summary(se)
       ID               X                Y               CA_m2      
 Min.   :  1.00   Min.   :676744   Min.   :5091658   Min.   : 6.48  
 1st Qu.: 30.25   1st Qu.:676754   1st Qu.:5091661   1st Qu.:34.12  
 Median : 59.50   Median :676776   Median :5091663   Median :43.74  
 Mean   : 59.50   Mean   :676782   Mean   :5091663   Mean   :45.66  
 3rd Qu.: 88.75   3rd Qu.:676802   3rd Qu.:5091666   3rd Qu.:56.17  
 Max.   :118.00   Max.   :676832   Max.   :5091669   Max.   :95.58  
# What is the product of the function? Is it a raster or vector?
se
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 118, 4  (geometries, attributes)
 extent      : 676744.4, 676837.1, 5091659, 5091747  (xmin, xmax, ymin, ymax)
 coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632) 
 names       :    ID         X         Y CA_m2
 type        : <int>     <num>     <num> <num>
 values      :     1 6.768e+05 5.092e+06 59.94
                   2 6.768e+05 5.092e+06 82.62
                   3 6.768e+05 5.092e+06 54.27
plot(se, axes = T)

### Let´s overlay the image and the product of segmentation (run both lines)
plot(forest)
plot(se, axes = F, add = T, main = "Lidar image of forest with segmented tree-crowns")

Task 11: Visualise the segmentation result in Leaflet

You can probably do all of this yourself, but here is a hint about projecting the SpatVector, just in case:

# What are we reprojecting and what to? geographical coordinates or?
se  # it is a Spatial object
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 118, 4  (geometries, attributes)
 extent      : 676744.4, 676837.1, 5091659, 5091747  (xmin, xmax, ymin, ymax)
 coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632) 
 names       :    ID         X         Y CA_m2
 type        : <int>     <num>     <num> <num>
 values      :     1 6.768e+05 5.092e+06 59.94
                   2 6.768e+05 5.092e+06 82.62
                   3 6.768e+05 5.092e+06 54.27
crs(se)  # what is its crs?
[1] "PROJCRS[\"WGS 84 / UTM zone 32N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 32N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],\n        BBOX[0,6,84,12]],\n    ID[\"EPSG\",32632]]"
# Project the SpatialVector using terra library
se4326 <- project(se, "EPSG:4326")

# ...then we can combine them with leaflet
library(leaflet)
leaflet() %>%
    addTiles() %>%
    addRasterImage(forest) %>%
    addPolygons(data = se4326, weight = 1, color = "black")  # Tadaa

Can you figure out where this forested landscape is from? And how many tree crowns have you just detected in how big an area? What is the average tree density?

Task 12: a more demanding segmentation example with a bigger LIDAR image!

And because we have it in the data/ folder, try with this larger example.

r <- rast("../data/myDem_subset.tif")
r
plot(r)

se_large <- itcIMG(r, epsg = 25829, ischm = T)  # can be a bit slow (2-3mins)
summary(se_large)
plot(r)
plot(se_large, axes = F, add = TRUE)

# Adjust 'th' argument for excessive capture of small growth
se_large5 <- itcIMG(r, epsg = 25829, th = 5, ischm = T)  # th - how low should algorithm be looking for canopy

# Plot raster with axes (plotting SpatVectors with terra tends to duplicate
# axes)
plot(r, axes = TRUE, main = "Forest LiDar with tree-crown overlay")

# Add vector without drawing axes again
plot(se_large5, border = "white", axes = FALSE, add = TRUE)

Again, what is the area of your raster and how many trees or bushes did you detect? Where do you think you are?

Want to see the result in Leaflet?

# Write a SpatVector to a Shapefile writeVector(se_large5,
# '../data/itcTrees_subset.shp', filetype = 'ESRI Shapefile')

# Load shapefile as SpatVector
rse <- vect("../data/itcTrees_subset/itcTrees_subset.shp")

# Reproject to WGS 84 (EPSG:4326)
rse4326 <- project(rse, "EPSG:4326")

# Control question: where is this landscape from?
r <- rast("../data/myDem_subset.tif")
leaflet() %>%
    addTiles() %>%
    addProviderTiles("Esri.WorldPhysical") %>%
    # addProviderTiles('Esri.WorldImagery') %>%
addRasterImage(r) %>%
    addPolygons(data = rse4326, weight = 1, color = "black")  # Neat :)

Task 13: Classify green areas and find the tree crowns in downtown Aarhus

a <- rast("../data/Aarhus_1m.TIF")
plotRGB(a)

library(RStoolbox)

# unsupervised classification with 3 classes
uc <- unsuperClass(a, nClasses = 4)

# plot result using ggr (ggplot for rasters)
library(ggplot2)
ggR(uc$map, geom_raster = T, forceCat = T) + scale_fill_manual(values = c("grey",
    "darkgreen", "sandybrown", "blue"))

Try supervised classification of Aarhus landuse with training data

library(sf)
crs(a)
[1] "PROJCRS[\"ETRS89 / UTM zone 32N\",\n    BASEGEOGCRS[\"ETRS89\",\n        DATUM[\"IRENET95\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1,\n                    ID[\"EPSG\",9001]]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]],\n    CONVERSION[\"Transverse Mercator\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"northing\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"
train <- read_sf("../data/train_aarhus.geojson") %>%
    st_transform(crs = crs(a))
st_crs(train)
Coordinate Reference System:
  User input: PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        DATUM["IRENET95",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]] 
  wkt:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        DATUM["IRENET95",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
unique(train$Type)
[1] "water"    "forest"   ""         "building" "paved"   
train <- train[-7, ]
train$Type <- factor(train$Type, levels = c("building", "forest", "paved", "water"))

# plot input data
ggRGB(a, r = 3, g = 2, b = 1, stretch = "lin") + geom_sf(data = train, aes(fill = Type)) +
    scale_fill_manual(values = c("yellow", "darkgreen", "sandybrown", "blue"))

# fit random forest (splitting training into 70\% training data, 30\%
# validation data)
sc <- superClass(a, trainData = train, responseCol = "Type", model = "rf", tuneLength = 1,
    trainPartition = 0.7)

# print model performance and confusion matrix
sc$modelFit
[[1]]
  TrainAccuracy TrainKappa method
1     0.8880996  0.8472588     rf

[[2]]
Cross-Validated (5 fold) Confusion Matrix 

(entries are average cell counts across resamples)
 
          Reference
Prediction building forest paved water
  building    667.8    5.6  66.4  23.6
  forest       22.6  326.6   9.6  13.2
  paved        83.2    2.6 482.4   3.0
  water        14.6    3.0   4.6 523.2
                            
 Accuracy (average) : 0.8881
# plotting: convert class IDs to class labels (factorize) and plot
r <- as.factor(sc$map)
levels(r) <- data.frame(ID = 1:4, class_supervised = levels(train$Type))
ggR(r, geom_raster = T, forceCat = T) + scale_fill_manual(values = c("sandybrown",
    "darkgreen", "yellow", "blue"))

Task 14: Optional: Spectral unmixing

Sometimes you want just a two-color/binary output, showing the breakdown on two main landuses (e.g. water vs urban, or forest vs paved, etc.) RStoolbox offers spectral unmixing by implementing the Multiple Endmember Spectral Mixture Analysis (MESMA) approach for estimating fractions of spectral classes, such as spectra of surfaces or materials, on a sub-pixel scale. You can adapt the workflow below for an image of your choice.

The following workflow shows a simple Spectral Mixture Analysis (SMA) with single endmembers per class, extracted from the lsat example image by cell id:

library(RStoolbox)
library(terra)

# to perform a SMA, use a single endmember per class, row by row:
em <- data.frame(lsat[c(5294, 47916)])
rownames(em) <- c("forest", "water")

# umix the lsat image
probs <- mesma(img = lsat, em = em)
plot(probs)

Interactively select raster IDs

library(terra)

# Plot the image (pick one band or RGB to visualize)
plotRGB(a, stretch = "lin")  # or RGB with plotRGB() if applicable

# Interactively click on the map (e.g., to pick endmembers)
clicks <- click(a, n = 4, cells = TRUE)  # 'n' = how many points to click

# View clicked cell IDs and coordinates
print(clicks)

The End

Similar approaches can be used when mapping socio-cultural phenomena in satellite imagery, from mosaicing images of night lights as proxies of economic performance, or detecting phenomena in the landscape such as burial mounds, growing urban sprawl, or tracing the outlines of scanned line drawings. (In the latter two you may need to base the classification on reflectance or edge detection rather than elevation.)